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We use the time dependent variational matrix product state (tVMPS) approach to investigate 
the dynamical properties of the single impurity Anderson model (SIAM). Under the Jordan- Wigner 
transformation, the SIAM is reformulated into two spin-1/2 XY chains with local magnetic fields 
along the z-axis. The chains are connected by the longitudinal Ising coupling at the end points. 
The ground state of the model is searched variationally within the space spanned by the matrix 
product state (MPS). The temporal Green's functions are calculated both by the imaginary and 
real time evolutions, from which the spectral information can be extracted. The possibility of using 
the tVMPS approach as an impurity solver for the dynamical mean field theory is also addressed. 
Finite temperature density operator is obtained by the ancilla approach. The results are compared 
to those from the Lanczos and the Hirsch-Fye quantum Monte- Carlo methods. 
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I. INTRODUCTION 

During the past twenty years, the dynamical mean 
field theory (DMFT)0^has been quickly developed into a 
powerful method to solve the strongly correlated models 
on the lattice (for a review see Georges et al.^^. DMFT 
maps the lattice models to the corresponding quantum 
impurity models subject to the self-consistency condi- 
tions. Unlike the normal static mean field approaches, 
DMFT keeps the full local dynamics induced by the lo- 
cal interactions. DMFT has been successfully applied 
to various correlation problems, such as the Mott tran- 
sition in the Hubbard modeP^ and the heavy fermion 
system^^. 

In DMFT one encounters the problem of how to effi- 
ciently solve the quantum impurity problems with self- 
consistently determined bath. The impurity solver can 
be regarded as the engine of DMFT, which influences 
the efficiency and accuracy of DMFT calculations. Since 
the invention of DMFT, many impurity solvers have 
been developed. With the development of the mod- 
ern computers, the essentially exact numerical methods 
have received much attentions. Among them the most 
used methods include the exact diagonalization (ED)'', 
Hirsch-Fye Quantum Monte Carlo (HFQMC^f^ and the 
numerical renormalization group (NRG)* . Most recently 
the continuous-time quantum Monte Carlo (CTQMC) 
solvei^^ has been introduced. These methods can be 
divided into two classes. ED and NRG^^ approaches are 
based on the Hamiltonains, while the QMC solvers are 
based on the Lagrangians (action). Generally the La- 
grangian approach is favored since the DMFT theories it- 
self is derived in the Lagrangian representation, and thus 
in the self-consistent process there is no need to map the 
continuous hybridization function to the discrete Hamil- 
tonians. However the Hamiltonian approaches have the 
merit that they work well at low temperatures and for 



real frequency which is more relevant to the experimental 
quantities, since most of the novel quantum phenomena 
in condensed matter physics happen at very low tem- 
perature. By the term of "impurity solver", one means 
not only to compute the ground state but also the whole 
spectral functions, which include the lower energy quasi- 
particle parts as well as the higher energy Hubbard bands 
in general. NRG approach could resolve exponentially 
small scales at the expense of the accuracy at intermedi- 
ate and high energies. 

The extensions of DMFT also call for the develop- 
ment of impurity solvers. In recent years, LDA+DMFT 
has been developed very quickly and successfully ap- 
plied to many system^, see Kotliar et al. and Heldl^ 
for the reviews of the recent developments and applica- 
tions. Since in the real materials there are usually more 
than one orbitals involved, therefore one needs to solve 
the quantum impurity problems with multi-orbitals effi- 
ciently. A second direction of the development is to study 
not the impurities but small clusters embedded in bath 
to capture the spatial fluctuations^^. Therefor one needs 
an efficient method to calculate the spectral properties 
of smaller clusters. The requirement of solving complex 
impurity problems (multi-orbital or cluster) puts hurdles 
on the ED and NRG approaches. 

A third direction of extension of the DMFT is to deal 
with the non-equilibrium system^^, which do not have 
the translational symmetry in the time domain. It re- 
quires the solver being able to work in the time domain, 
while to our knowledge most of the solvers work in the 
frequency (imaginary or real) domain. A recent attempt 
is made in the CTQMC approach^^. While this approach 
can handle arbitrary interaction strengths, it suffers from 
a dynamical sign problem which becomes severe at long 
time or for large bandwidth. 

Therefore it is urgent to develop an impurity solver 
working at zero temperature, which satisfies the following 
criteria, i) It can capture both the low energy quasi- 
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particle physics and the high energy Hubbard bands, ii) 
It is easy to be generahzed to multi-orbital or cluster 
cases, iii) It works with real frequency and gives the real 
time dynamical properties directly. 

In this paper, we develop an impurity solver satisfy- 
ing all these criteria based on the time-dependent vari- 
ational matrix product states (tVMPS) approach. The 
variational approach directly attacks the strongly corre- 
lated problems by an educated guess of the many body 
wavefuncitions. By introducing more variational pa- 
rameters one could enlarge the dimensions of the vari- 
ational space but it makes the variation process harder 
in general. For example, in the case of Gut z wilier varia- 
tional wave function, the evaluation of the expectation 
value is often done approximately by introducing the 
Gutzwiller approximation^^ Matrix product states 
(MPS) are states where the coefficients of the wave func- 
tion are a product of matrices depending on the local 
lattice site states. They are generated naturally from the 
NRG and DMRG calculations. Actually the latter ap- 
proach can be reformulated into the variational approach 
within the space spanned by the MPS^^. Time evolution 
algorithrrP^ based on MPS was first proposed from the 
quantum information perspective, and then been trans- 
lated into the language more access to the many-body 
theorists^^*^, see Verstraete et al.l^, Garcia- Ripoll ^-^ 
and Schollwoeck and White ^"^ for reviews of the time 
evolution algorithm of the MPS. 

There are some previous attempts using the VMPS ap- 
proach to study the quantum impurity problems. Weich- 
selbaum et al.^^ studied the spectral property of SIAM 
in the presence of a magnetic field using the correc- 
tion vector approach, which is first proposed in the con- 
text of dynamical density matrix renormalization group 
(DDMRG)^^ Saberi et al.'^ performed detailed compar- 
ison of the VMPS and the NRG approach on the SIAM. 
Holzner et al. used the VMPS formalism to study the 
two-lead, multi-level Anderson impurity model. Along 
the line of developing fast impurity solvers, several of the 
authors have developed the Gutzwiller based impurity 
solvei^, which is suited for combining with the LDA 
and studying the low temperature properties of multi- 
orbital models. There were previous attempt of us ing 
the DMRG as the impurity solver of the DMFT.I^^ 

The organization of this paper is as follows. In Sec. 
[TT|we give an overview of the computation of the ground 
state of the SIAM in the VMPS formalism, where an 
unfold^ technique greatly reduces the computational ef- 
fort. In Sec. |III| we describe the algorithm for time evolu- 
tion in both real and imaginary times, by which one could 
calculate the Green's functions in the time domain. In 
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Sec|IV|we describe the fit and extrapolate scheme to ex- 
tract the spectral function from the real time data. In 
Sec. [V|we address the DMFT self-consistent loop where a 
fitting procedure from the continuous hybridization func- 
tion to the SIAM chain Hamiltonian is needed. Finally, 
we conclude the paper by making remarks on the pros 
and cons of the present solver and make an outlook for 
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Figure 1: (Color online) Unfold the SIAM chain. The SIAM 
chain is separated into two parts with different spins, they are 
connected at the leftmost end where the bold bond denotes 
the Hubbard interaction U. 



future developments. In the appendix we generalize the 
solver to finite temperature, where ancillary sites which 
play the role of the heat bath are introduced. 



II. VMPS APPROACH TO THE GROUND 
STATE OF THE SIAM 

The action of the single impurity Anderson model 
(SIAM) is 

SsiAM = r^ir rdrH(r)[(S,-/i)J,,,,+A(r,rO]c,(rO 
Jo Jo 

drlln^ni (1) 

/o 

It is a zero dimensional problem, and there only exists 
the fluctuation in the temporal axis, which is captured 
by the time dependent hybridization function A(r, r'). 
One could de-integral the action by introducing the non- 
interacting bath degree of freedoms. The electron hops 
into the bath, travels in it for a while and then hops back 
to the impurity site thus brings in the temporal nonlocal 
correlations. The de-integral process is not unique, and 
in the DMFT context, both the star and the chain ge- 
ometry of the bath degree of freedom have been studied. 
In the present study we reformulate the SIAM into a one 
dimensional chain with the nearest neighbor hopping, see 
Figure [T] In the framework of DMFT the Hamiltonian 
parameters of the chain are determined self-consistently, 
and generally has no translational invariance. We deal 
with the chain with typical chain length of 10 to 20 sites 
in the present study. 

By introducing the Jordan- Wigner transformation 
(JWT) for the spin up and down fermions separately, 
^ii^r) = {Ylk<iPka) Cia where pka = (-l)""^- is the 
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Figure 2: (Color online) (a). The A^^^'^^ has three indices, ak 
is the physical index run from 1 to d = 2 for SI AM. i (j) is 
the visual index runs from 1 to Xfc (Xfc+i)- They control the 
precision of the VMPS calculation, (b). An operator acts on 
a single site (c). A MPS where all of the connected bonds 
have been contracted. 



JW sign, one could unfold the SIAM chain into two spin- 
^ XY chains. They are coupled at the end point by the 
longitude Ising coupling due to the Comlomb interaction. 
The onsite energy of each site is mapped into the local 
magnetic filed along the z-direction. Since the hopping 
is between the nearest neighbors, the JWT sign does not 
shows up in the calculation of the ground state. A simple 
division of the Hamiltonian into terms acting on the odd 
and even bonds are possible, see Sec |III[ 
In terms of 1/2 spins, the Hamiltonian is 
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HsiAM = Himp + Hi)ath + Hhyh 

Nbath 
Nbath-l 

i=l 

Hhyh = Ji{s^s~ ^r^r^ ^h.c.) 

(2) 

where sf = ^""^^^^ are the ladder operators for spin 
one half, s and r denote the spin up/ down fermion 
operators. 

The matrix product state (MPS) H^mps) = 

E.,..,......„1V(4'^U^l...Al^-l)|<Ti,<T2,...,<T^). a, = 

1, is the physical index, where d = 2 is the dimen- 
sion of the local Hilbert space, is a matrix with 
the visual dimension Xk x X/e+i- There is an efficient 
scheme based on the transfer matrix to handle them in 
low dimensions, i.e. to calculate the overlap between two 
MPSs, and the expectation value of an operator over two 
MPSs, the operation of a local operator on a given MPS. 
They are shown schematically in Fig |3] and see Verstraete 
et al.— for a reference. 

We minimize the expectation value 



Figure 3: (Color online) (a). The expectation value of the 
operator H over the MPS. (b). Tracing out all the indices on 
the sites other than k gives the effective Hamiltonian Heff 
on the k-th site. (c). The minimization problem on site k is 
equivalent to the eigen-problem of the effective Hamiltonian 



{i^MPs\HsiAM\i^MPs) within the space spanned 
by the normalized MPS. Such a problem is solved very 
efficiently by the alternating least squares scheme, in 
which we perform the optimization of the matrices site 
by site. In each step, we fix all but the matrix on the 
current site, see Figure [3] By contract all of these indices 
we define the effective Hamiltonian on the current site. 
The minimization problem becomes quadratic and is 
equivalent to an eigen-problem of the size Xk x X/e+i x d. 
Since only the ground state is needed, we solved it by 
the Lanczos method. The key step in the Lanczos step 
is the computation of the matrix-vector multiplication: 

Where H = can be parallelized over each term 

in the Hamiltonian. 



III. TIME EVOLUTION 

To calculate the retarded Green's function G^{t) = 
—i6{t){GS\{C(j{t),c'l}\GS), one first applies the cj. or c^r 
on the ground state from the previous VMPS calcula- 
tions: |(^(0)) = cl\GS), |x(0)) = c^\GS). Then let states 
evolute in real time to get = e-'^^^-^^^^cllGS) and 

\x{t)) = e'^^-^^^^c^\GS). Thus the retarded Green's 
function can be calculated as 

Giit) = -iemmim + ixmxm (3) 

The operation of the local operators on the ground 
state MPS is 

l<k myk 

(4) 

Where a string of JW signs and di dx d matrix 5+ act on 
the k-th site. After the multiplication, the new state is 
still represented as the the MPS of the same rank. 

The real-time evolution is performed by first split- 
ting the evolution operator into small pieces e~*^^ = 
^^-iHAt^M^ where At = ^, then applying the Trotter 
decomposition. To the second order one has e~*^^^ = 

^-iHoAt/2^-iH,At/2^-iHz At ^-iHe At/2 ^-iHo At/2 ^ 

0((At)3). Where H = Ho + ^ Hz- and Ho 
are two particle terms act on even and odd bonds and 
Hz includes the single particle terms. Since each term 
within Hoc^ (o^ = e,o, Z) commutes one could apply the 
evolution operator Ua{^t) = e~*^"^^ associate with 
each bond one after another without leading to further 
errors. In the particle-hole symmetric case, the onsite 
energies of the chain are zero and the chemical potential 
li = U/2 can be put into the interaction term. Thus the 
single particle evolution operator Uz = I 

The operation of the single site evolution operator Uz 
is similar to the applying of the creation/annihilation op- 
erators, Eqn |4] The operation of the two site evolution 
operators Ue or Uo is shown schematically in Figj4] The 
physical indices for the two sites are first merged and ex- 
changed, then the d'^ x d'^ operator is separated into two 
matrices Ui and U2 by the singular value decomposition 
(SVD). Each of them only carries the physical indices of 
one site. In both cases, the evolution operators are writ- 
ten as the matrix product operator (MPO)^^, which has 
four indices, two for the visual dimensions and two other 
for the physical indices, see Figj2] Generally after apply- 
ing them onto the target states, the resulting states could 
not be written as the MPS with same visual dimension as 
before. Thus one needs to project the state into the orig- 
inal MPS space. It is done variationally by minimizing 
the norm: A/" = — Ua{^t)\(j){t))\\^ and approximately 
one has \(j){t + At)) — The minimization procedure 
can be carried out following the similar procedures as 
in Sec|n[ It is performed every step after the applying 
of the evolution operator on the original state. In the 
course of evolution, we increase the bond dimension to 
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^^igure 4: (Color online) Writing the two particle evolution 
operator into a MPO. (a). The physical indices of the four 
indices tensor are merged, giveing a x matrix, (b). 
The order of the indices are exchanged, (c). SVD of the 
matrix separates it into two parts (d). Rearrangement of the 
physical indices gives the MPO representation of the two-sites 
evolution operator. 



maintain the smallest eigenvalue of the reduced density 
matrix larger than a threshold Wmin — 10~^^. 

See Figure [5] for the real time Green's function for dif- 
ferent interacting strength. The retarded Green's func- 
tion G^(t) evolutes from the Bessel function at the non- 
interacting limit to the cosine function at the atomic 
limit. In between one has the oscillated decaying curves, 
and the period of the oscillation gives the frequency of 
the Hubbard band. The normalization of the density 
of states (DOS) is automatically fulfilled since it implies 
that G^(t = 0) = —1. There may be a concern on the 
orthogonality catastrophe (OC) F^^ which states that in 
the thermodynamic limit local perturbation leads to com- 
plete reconstruction of the ground state of a fermionic 
system in such a way that the overlap of the "old" and 
"new" ground-state wave functions is proportional to 
. But si nce w e are dealing with finite N here, the OC 
is irrelevant j ^^l^^ i By comparing the real time Green's 
function with different bath sites, we notice that even 
with small number of sites, one could reproduce the re- 
sulting spectral information for the thermal-dynamical 
limit, as long as the wavefront created by the adding (re- 
moving) an electron on the impurity site does not reach 
the boundary. 

Imaginary time Green's function could also be calcu- 
lated similarly, in which we perform the evolution algo- 
rithm along the imaginary time axis. Imaginary time evo- 
lution technique could also be used to search the ground 
state. The imaginary time Green's function G(r) is a de- 
cayed function from n^r — 1, and it is decaying faster for 
larger interacting strength U . For the particle- hole sym- 
metric case, the decaying form approaches single mode 
decay —\e~^^ in the atomic limit. 
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Figure 5: (Color online) The zero temperature retarded 
Green's function G"^(t) for a 20-site SIAM chain. Where 
Sn = 0,7n = 0.5, the interacting strength U — 0,1,2,3,4,5 
from the bottom to top. The diamond shaped dots indicate 
the analytical result for U — 0. Finite size effect of the chain 
does not show up for this time scale, so the Fourier transform 
of the real time Green's function gives spectral functions. 
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amount of raw data. Generally there are two approaches 
to avoid the overboadening: the linear predic tion (LP) 
and the fitting and extrapolate (FE) procedur^^^EHl. We 
leave the LP approach for further investigation and adopt 
the FE approach which we believe is more under control. 

We fit parts of the raw data with Acos{uj{t — tQ))/t'^ + 
Be~^^^ where A, B^to, f3 are fitting parameters. Then 
the function is extrapolate to very large time. By this 
method we have taken use of most of the raw data. A 
comparison of DOS from directly FT and after "fit and 
extrapolate" is shown in Fig|6] Even from directly FT 
data one could recognize zero component Kondo peak 
and the development of Hubbard band as U increases. 
However, the Fridel sum rule which states the pinning of 
the A{0) = is violated. Actually, the conservation 
of zero frequency peak relies on the conservation of the 
areas enclosed by the G^{t) curve. This means the res- 
olution of the Kondo peak relies on very accurate long 
time behavior of the Green's function. By the FE proce- 
dure before FT, this conservation could be fulfilled, Figj6] 
However, the fitting ansatz inevitably introduces prior 
knowledge about the spectral and gives much sharper 
peaks in present case. Note that causality is violated 
(negative dos) in /7 = 1 case. 

We find it is more stable to extract the spectral 
function from the G(r) by maximum entropy method 
(MEM). Since the data is free from statistical error, the 
procedure is more stable and less involved than the con- 
tinuation of Monte Carlo results. As stated in early work 
of Naef et ai"^^ , although the continuation approach may 
be still insufficient for resolving specific line shapes, it 
is appropriate for identifying gaps and peak positions, 
which is more relevant in DMFT studies. So in present 
paper MEM is still adopted to get DMFT converged spec- 
tral functions. 



V. DMFT SELF-CONSISTENT 



Figure 6: (Color online) DOS of a SIAM for U = 1,2,3,4 
from directly FT (full line) of real time data (FigjsJ and fit 
the long time behavior and then perform FT (dash line). 



IV. EXTRACT SPECTRAL FUNCTION 

The zero temperature spectral information can be ex- 
tracted from the imaginary time Green's function G(t), 
by the maximum entropy continuation similar to the pro- 
cedures in the PQMC approach to the SIAM^^. Or one 
could get the spectral information from the Fourier trans- 
form of real time Green's function G^{t).^^ 

In the Fourier transformation (FT) approach, to avoid 
the negative DOS, generally one multiplies the time do- 
main data with a window function {W{t) = e~^^^rnax ^ is 
used here) which is a decayed function from zero time to 
the large time cutoff, but this method inevitably broaden 
the peak in the spectral function and also drops a large 



After getting the temporal Green's function, one has 
essentially solved the impurity problem. One should plug 
it into the DMFT self-consistent loop to produce the 
SIAM chain Hamiltonian for the next iteration step. For 
demonstration purpose we consider the case of the Bethe 
lattice where the self-consistent procedure is greatly sim- 
plified. By multiplying ^ with G{t) one gets the hy- 
bridization functions A(r), where D is the half-bandwith 
of the semi-ellipse DOS. For the general lattice, the self- 
consistent loop should be 1) Using Fourier transform or 
MEM to get the Green's function in the frequency do- 
main, and do the self-consistent following the conven- 
tional routine. This kind of self-consistent procedure was 
adopted in PQMC solvei^^. or 2) Just do all of the com- 
putation in the time domain. In that case, one needs to 
do inverse of the Green's function matrix in time domain 
and convolution of non-interacting DOS with each matrix 
element. It is a formidable task but still could be done, 
as one encountered in the non-equilibrium DMFT.^^- 
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Figure 7: Flow diagram of DMFT self-consistent calculations. 
tVMPS solves Green's function for a SIAM chain Hamilto- 
nian in temporal domain. Self-consistent process gives new 
hybridization function A(r), from which the Hamiltonian pa- 
rameters {£i,7i} are fitted, thus closes the loop. (a). For 
general lattice, one performs Fourier transform or analytical 
continuation to translate G{t) to frequency domain, and then 
determines new hybridization function, information of non- 
interaction DOS enters in this step. (b). The self-consistent 
loop simplifies for infinite dimension Bethe lattice, where Mat- 
subara Green's function G{r) directly determines the new hy- 
bridization function. Spectral function A{uj) could be gotten 
by MEM continuation of the converged Matsubara Green's 
function G{r) or FT of retarded Green's function G^{t), in 
the latter case, real time evolution using converged Hamilto- 
nian parameters is firstly performed. 



The key problem then is how to determine the Hamil- 
tonian parameters of the SIAM chain from the contin- 
uous hybridization function A(r). Similar problem is 
encountered in other Hamiltonian based solvers such as 
ED, NRG and DDMRG approaches. In the ED ap- 
proach the step is determined by the conjugated gradi- 
ent minimization^ or continued fraction expansion^^ of 
the Green's function in the frequency domain, but due 
to the finite size of the effective bath that can be dealt 
with, one has to truncate and leads to errors in this step. 
In the NRG approach the logarithmic discretization pro- 
cedure is adopted^ while in the DDMRG approach, 
no logarithmic energy separation is assumed, so a direct 
tri-diagonalization scheme can be used, i.e. first fit the 
hybridization function with a star shaped Hamiltonian, 
then using Lanczos method to tri-diagonalize the Hamil- 
tonian, and the diagonal/off diagonal matrix elements are 
the onsite/hopping parameters for a chain Hamiltonian. 

In the present case, we fit A(r) with the Hamiltonian 



(5) 



A(r) is the hybridization function for a noninteracting 
SIAM chain with onsite energy Si and hopping amplitude 
7„ i = 1,2,. ..math- A(r) = Eni\Uii\^e-^'^[0{r){ni - 
1) + 0{—r)ni], where U and E are eigenstate and eigen- 
value of the noninteracting Hamiltonian Hbath , satisfying 
HbathU = UE. rii is the occupation number on level /, 
and ni = 0{—Ei) for zero temperature. We use conjugate 
gradient method to search within the parameter space 
spanned by Si and ji which minimize Xa- Generally 
there are 2Ni)ath parameters, while in the particle- hole 
symmetric case all of the Sn is kept as zero and there- 
fore one has Ni^ath fitting parameters. Since we are deal- 
ing with relatively larger number of bath sites comparing 
with ED approach and the function to be fitted is more 
well behaviored (different from the fitting G{iujn) in fre- 
quency domain, one has A(r) which is more smooth and 
featureless), the fitting is more reliable. In the present 
study, the average fitting error can be reduced to 10~^ 
per data point. 

After the DMFT loop converges, one could get various 
information from the converged Green's function G{r) 
and the MPS. e. g. the kinetic energy and the double 
occupancy rate. See Figj9] for the comparison of the re- 
sults from the tVMPS and the Lanczos solver at zero 
temperature. The long time asymptotic behavior for the 
interaction strength U larger than the critical Uc shows 
the presence of the charge excitation gap. One could also 
recognize the metal-insulator transition clearly from the 
DOS on the Fermi surface ^4(0) = — lim/3^00 f ^( 2 ) 
the suppression of the double occupancy. The real fre- 
quency information can be extracted from the analytical 
continuation approaches, such as the Pade approximation 
or the MEM. The MEM inverted the integral transfor- 
mation —G{t) = j'^^K{uj,T)A{uj), where symmetrical 

kernel K{uj,r) = ^- — ^J^t''^-^^^^ is used to persevere the 
particle hole symmetry A{uj) = A{—uj), /3 is a fictitious 
temperature /3 = 128 in our calculation. The results are 
shown in FigfTo) 

It is also possible to do the DMFT self-consistent in the 
real time domain, but in the DMFT self-consistent loop 
we adopt the imaginary time Green's function which is 
more stable. And since all of the evolution operators are 
real, this choice would reduce the numerical effort. Once 
convergence is achieved, a real time evolution could be 
redone to get the retarded Green's function G'^{t) and 
the spectral function could be extracted from it with ap- 
propriate FT method. 



VI. SUMMARY AND OUTLOOKS 

Here we discuss various sources of error and computa- 
tional effort of the tVMPS solver. First due to the cut- 
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Figure 8: (Color online) Main figure: The DMFT con- 
verged zero temperature Green's functions by the tVMPS 
solver with 18 sites. Where the interaction strength U = 
0, 0.8D, 1.6D, 2 AD, 3. 2D, A.^D from bottom to the top. The 
long time asymptotic behavior of G{t) differentiates the 
metallic and the Mott insulator states. In set: The density of 
states on the Fermi surface A(0), which are compared with the 
results from a Lanczos solver with 8-site. Slightly dropping 
of A(0) in the metallic side is due to finite /3, which prohibits 
to resolve very small frequency. 




U/D 



Figure 9: (Color online) Physical quantities at T = from 
the tVMPS solver compared with the Lanczos solver. Up- 
per panel: The average kinetic energy per site, Lower panel: 
The average double occupancy per site versus the interacting 
strength U/D. 



off of the matrix size in the MPS representation, there 
are errors in the ground state and the evolution proce- 
dure. Also there are Trotter decomposition errors from 
the discretization of the evolution operators. Finally, er- 
rors also come from the mapping from the hybridization 
function to the chain Hamiltonian. The latter two errors 
can be easily reduced by adopting higher order Trotter 
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Figure 10: (Color online) The DOS extracted from the con- 
verged G{t) by MEM. 



decomposition and by increasing the number of the fit- 
ting parameters, i.e. the length of the chain. However, 
the first error is more severe and prohibit long time evo- 
lution. It is known that the algorithm will breakdown 
at runaway time, which approximately shows the loga- 
rithmic dependence on x- At the runaway time the en- 
tanglement of the states increases to a value which could 
not be sufficiently represented by the given finite trunca- 
tion dimension. The computational effort of the tVMPS 
solver is mainly determined by the truncation dimension 
X and is insensitive to the Hamiltonian parameters of the 
problem. 

We then summarize the limitations due to present im- 
plementation of the algorithm, and point out possible 
solutions. First, the algorithm does not have a good res- 
olution on the Kondo resonant peak, i.e. it does not pin 
down at the noninteracting value as U increases. This 
could be cured by u sing the logarithmic discretization 
in the bath, an d has attained some success in using 
this technique to resolve Kondo effect in the transport 
properties. Second, to reduce computational effort we 
limit chain length only slightly larger than that acces- 
sible to exact diagonalization, but by using conserved 
quantity such as total particle number, magnetization or 
even non-abelian SU (2) symmetry of SIAM, much larger 
system size is accessible within the same computational 
time. Third, accurate determination of low temperature 
Green's function needs long time evolution without hit- 
ting the runaway time. Whether the tVMPS could pro- 
duce result that competes to the Monte Carlo method 
should be answered by future investigation. 

To conclude, we have shown that the tVMPS approach 
could work at zero temperature as well as at finite tem- 
peratures. By the post-selfconsistent real time evolution, 
one could also get the real time dynamics, from which the 
real frequency can be extracted. These merits make it a 
promising solver for DMFT to investigate the low tem- 
perature properties of the strongly correlated systems. 
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As an outlook, the time- dependent VMPS solver re- 
ported here can be generalized to complex impurities, 
for example, the multi-orbital impurity problems, and 
the quantum cluster problems. Different from the ED 
based solvers, the generalization of the present solver to 
multi-band or clusters case does not encounter the ex- 
ponential increasing of states. Therefore tVMPS is a 
promising solver to work within LDA+DMFT^^^ to in- 
vestigate the realistic materials and a cluster solver works 
within the cluster-DMFT formalisms."^. The Trotter er- 
rors can be avoided even by adopting different evolution 
algorithm, i.e., the Lanczos dynamics^^ and the time step 
targeting^^. Different contract schemes (transverse con- 
tract) in the time evolution are reportecP^. New algo- 
rithm for the finite temperature based on the minimally 
entangled typical quantum states is reported in^^ and 
can be used to replace the ancilla approach. The present 
solver works in real time domain, and thus may be use- 
ful in the non-equilibrium DMFT^^. In an independent 
direction, the impurity solver developed here could also 
be used to study the time-depend ent p henomena in the 
quantum transport in nanode vices. ^^^^ 
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Appendix A: Finite Temperature Algorithm 

In this section, we discuss the finite temperature algo- 
rithm based on the ancilla approach. The ancilla 
approach replaces the mixed state needed by comput- 
ing the thermodynamical averages by the pure state of 
an enlarged system. The enlarged system is constructed 
simply by the original physical system 1-L and the identi- 
cal copy of it, whose Hilbert space is denoted by A. For 
the present case, the ancillary sites are added parallel 
to the SIAM chain, make the system resemble a ladder, 
and the Hamiltonian acts only on the physical sites, see 
FigjlT] One first prepares the states. 



Tr^ I -0/3) (-0/3 1 and the partition function Z = {tp isltp 13) . 
The finite temperature correlation function can be eval- 




(Al) 



where the first (second) cr^ G 1-L{A) denotes the state 
of the i-th physical (ancillary) site. By the purifica- 
tion process one gets {ipfs) = e~^^^'^\ipQ) ^ from which one 
could get the finite temperature density operator pf3 = 



uated as Tr-^ (p^c(r)c''' 



1^/^/3). So once 



one gets IV^/^), one could apply the evolution algorithm 

Physical 



Ancilla 



Figure 11: (Color online) Introduce the ancillary sites parallel 
to the SIAM chian. The ancillary sites act as the heat bath. 
When been traced out they give the thermodynamic averages. 
The presence of the ancillary site enlarges the local Hilbert 
space. The Hamlitonian only acts on the physical chain. 



similar to the Sec III to compute the Green's function in 
the imaginary as well as real times. 

The imaginary time Green's function G(r) with r = 
to j3 is computed from the evolution algorithm. For 
the particle-hole symmetrical case, one only needs half 
of the data (r < /3/2) due to the symmetric property 
G{t) = G{j3—T). For general cases, we evolute G{t) from 
to ±13/2 and then use G(/3/2 < r < /3) = -G(-/3/2 < 
r — /3 < 0) to restore the imaginary time Green's function 
from to /3. This approach reduces the computational 
efforts and also the accumulated errors in the long time 
evolution. (The ending point is conserved since in the 
present approach it is equal to (?/^o|ce~^^/^c'''e~^^/^|'?/^o) 
) To compute the Green's function at low but finite 
temperatures, one needs to perform long time evolution, 
which is not stable due to the accumulation of the Trot- 
ter error and the truncation of the Hilbert space. The 
error is more severe in the finite temperature case, be- 
cause of the presence of the exponent growth factor e^^ . 
This hinders the investigation of very low temperatures 
but can be circumvented by choosing alternative evolu- 
tion algorithmP^lEHl. These possibilities demands further 



investigating. See Fig 12 for a comparison of the finite-T 
Green's function for a SIAM chain from the tVMPS and 
the HF-QMC approach. 

There is an alternative way which does not group the 
physical and the ancillary sites into the supersite. This 
way reduces the size of the local Hilbert space, but the 
non next-nearest interaction will hinder the simple Trot- 
ter decomposition based time evolution algorithm. After 
all, the time evolution can be done by the Lanczos dy- 
namic or the time step targeting which does not assume 
the next-nearest interactions. 
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Figure 12: (Color online) The /3 = 16 imaginary time Green's 
function for a 4-site SIAM chain calculated by the tVMPS 
approach (lines) compared with the HF-QMC results (hollow 
dots). En = 0, 7n = 0.5, U = 0, 1, 2, 3, 4, 5 from bottom to top. 
r > /3/2 data are calculated from the — /3/2 < r < data by 
G(t) — —G{t — P) to avoid long time evolutions. 
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